Giant anharmonicity and non-linear electron-phonon coupling in MgB 2 ; 
Combined first-principles calculations and neutron scattering study 
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We report first-principles calculations of the electronic band structure and lattice dynamics for the 
new superconductor MgB2 . The excellent agreement between theory and our inelastic neutron scat- 
tering measurements of the phonon density of states gives confidence that the calculations provide 
a sound description of the physical properties of the system. The numerical results reveal that the 
in-plane boron phonons (with E2 9 symmetry) near the zone-center are very anharmonic, and are 
strongly coupled to the partially occupied planar B a bands near the Fermi level. This giant anhar- 
monicity and non-linear electron-phonon coupling is key to quantitatively explaining the observed 
high T c and boron isotope effect in MgB2 



PACS numbers: 63.20.Ry, 63.20.Kr, 61.12.-q, 74.25. Jb, 74.25.Kc 



The recent discovery of superconductivity at 40 K in 
the MgB 2 binary alloy system M has triggered enormous 
interest in the structural and electronic properties of this 
class of materials. The system has a very simple crys- 
tal structure H, where the boron atoms form graphite- 
like sheets separated by hexagonal layers of Mg atoms 
(see inset to Fig. 1). Our pscudopotential plane wave 
band structure calculations show that the bands near 
the Fermi level arise mainly from the p XjV a bonding 
orbitals of boron, while the Mg does not contribute ap- 
preciably to the conductivity, in good agreement with 
initial reports from other groups |^-^|- In the case 
of graphite these a bands are full, but for MgB2 they 
are partially unoccupied, creating a hole-type f| con- 
duction band like the high-T c cuprates. In contrast 
to the cuprates, however, the normal-state conductivity 
is three-dimensional in nature instead of being highly 
anisotropic, thus eliminating the "weak-link" problem 
that has plagued widespread commercialization of the 
cuprates. The normal-state conductivity is also one 
to two orders-of-magnitude higher than either the Nb- 
based alloys or Bi-based cuprates used in present day 
wires, and this feature combined with low cost and easy 
fabrication |^,^| could make this class of materials quite 
attractive for applications. 

From a fundamental point of view the central question 
is whether the high T c in this new system can be under- 
stood within the framework of a conventional electron- 
phonon mechanism, or a more exotic mechanism is re- 
sponsible for the superconducting pairing. The observed 
boron isotope effect jl0| argues for an electron-phonon 
mechanism, while the positive Hall coefficient || sug- 



lations of the lattice dynamics (and electronic) calcula- 
tions for MgB 2 . Excellent agreement is found between 
theory and experiment. More importantly, the numeri- 
cal results demonstrate that the in-plane boron phonons 
near the zone-center (with E2 S symmetry at T) are very 
anharmonic and strongly coupled to the partially occu- 
pied conduction band (planar B a bands) near the Fermi 
level, providing the large electron-phonon interaction in 
this system. Due to these strong anharmonic phonon 
modes, the electron-phonon coupling is non-linear, pro- 
viding the essential ingredient to explain the high T c and 
boron isotope effect in MgB2. 

The results of the neutron measurements and calcu- 
lations for the phonon density of states are summarized 
in Fig. 1. For the experimental data, the 5g polycrys- 
talline sample was prepared in the usual way with the n B 
isotope to avoid neutron absorption problems |T^ ], and 
was characterized by magnetization, x-ray, and neutron 
diffraction. Inelastic neutron scattering measurements 
were made to determine the generalized phonon density 
of states (GPDOS), which is the phonon density of states 
weighted by the cross section divided by the mass of each 
atom. The data were collected from 7 K to 325 K on the 
Filter Analyzer Spectrometer in the range 5-130 meV, 
and on the Fermi chopper instrument for energies of 0.5- 
30 meV. Additional details of the experiment and anal- 
ysis can be found elsewhere |13 3. The experimental 



gests similarities with the cuprates 11 1. To answer this 
question, we have carried out inelastic neutron scatter- 
ing measurements of the phonon density of states, and 
compare these results with detailed first-principles calcu- 



data indicate two bands of phonons, one below 40 meV 
corresponding primarily to the acoustic phonon modes, 
and one above 50 meV that mostly involves the boron 
motions. The high-energy portion of the spectrum shows 
a sharp cutoff at about 100 meV. The data are in basic 
agreement with the results of Osborn et al. jl5| , although 
the present data were obtained with better overall energy 
resolution. Our results do not agree with Sato et al. [jl6| , 
who reported a feature at 17 meV in the GPDOS. We did 
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not observe this feature at any temperature and conclude 
that it is not intrinsic to the GPDOS. We also have mea- 
sured the temperature dependence of the spectrum and 
find no substantial changes in the basic features from 7 
K up to 200 K, while a very modest softening of some of 
the modes was observed in going to 325 K. 
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Fig. 1 (a) Generalized phonon density of states for MgB2 
(top curve) determined from inelastic neutron scattering mea- 
surements, and the calculated phonon density of states (DOS) 
(bottom curve). The inset shows the unit cell of MgB2 along 
with the Brillouin zone and the high-symmetry directions, 
(b) Calculated dispersion curves of all phonons and the cor- 
responding DOS (right panel). At the zone center (F), the 
symmetries of the modes are also indicated. The modes indi- 
cated by the red-line dominate the electron-phonon coupling. 

In order to understand the origin of the features ob- 
served in the neutron data, we have carried out detailed 
first-principles calculations of the electronic band struc- 
ture and lattice dynamics. The calculations were per- 
formed using the pseudopotential plane wave method 
JT7| , within the generalized gradient approximation [ p"8| . 
We employed plane waves with an energy cutoff of 500 
eV, and the ultra-soft pseudopotentials for Mg and B 
p9[ . The total energy and forces converged within 0.5 
meV/atom and 0.01 eV/A, respectively. Brillouin-zone 



integrations |^(| were carried out using dk = 0.02 A -1 , 
generating 15 x 15 x 11 and 9x9x6 k-points for 1 x 1 x 1 
and 2x2x2 supercells, respectively. The lattice param- 
eters a and c were relaxed before the lattice dynamics 
calculations, and the force constant matrix was obtained 
by the direct-force method using periodically repeated 
supercells pj|. The optimized lattice parameters are in 
good agreement with our experimental values within 1%. 
The calculated density of states obtained from these cal- 
culations is shown in Fig. la, and we see that the agree- 
ment between experiment and theory for the energies of 
the modes is excellent. The agreement for the intensi- 
ties is also good considering that the observed spectrum 
is a weighted density of states rather than the actual 
DOS, and the wave vector averaging may include some 
coherency effects [^3|. The calculated phonon disper- 
sion curves along the high-symmetry directions of the 
Brillouin-zone (see inset to Fig. 1) are shown in Fig. lb. 
All the modes are found to be quite dispersive as expected 
for the small size of the unit cell, but divide nearly com- 
pletely into acoustic (below « 40 meV) and optical (50- 
100 meV) bands, with the only exception being along the 
A-L direction. All the optical modes are weakly disper- 
sive along the T-A direction, reflecting the layered nature 
of the MgB2 crystal structure. 

There are four distinct phonon modes at the zone cen- 
ter r, as indicated in Fig. lb. The A2« and Bi g singly- 
degenerate modes involve only vibrations along the c 
axis; for Bi s the boron atoms move in opposite direc- 
tions while the Mg ion is stationary, while for the A.2 U 
mode the Mg and B planes move in opposite directions 
along the z-axis. The other two modes are doubly de- 
generate and involve only in-plane motions (along the x 
or y axes). For the Ei„ mode the Mg and B planes vi- 
brate in opposite directions along the x (or y) axis. The 
calculated energies of all zone-center phonons are listed 
in Fig. 2, and for the three types of modes just discussed 
there is good agreement between calculations reported by 
other groups For the E2 S mode, however, there are 

large discrepancies between the results reported for var- 
ious calculations. The boron ions for this mode vibrate 
in opposite directions along the x (or y) axis, with the 
Mg ions stationary. The discrepancies noted above are 
resolved when the strong anharmonicity associated with 
the in-plane displacements of the boron atoms is taken 
into account. 

To further investigate this behavior, we plot in Fig. 2 
the energy as the atoms are displaced by an amount u ac- 
cording to each of the four eigenmodes. For the Bi ff , Ai tt , 
and Ei„ modes the distortion energy can be very well ap- 
proximated by the harmonic expression E(u) = A2U 2 up 
to u/a = 0.065. It is unusual that the phonons remain 
harmonic at these large distortions, but this is consistent 
with the fact that we did not observe any substantial 
changes in the measured GPDOS with temperature. 
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Fig. 2 Energy curves as a function of boron displace- 
ment for each zone center frozen-phonon mode. The top 
panel shows that the Bi s , A2u, and Ei« modes are harmonic 
up to u/a ~ 0.065. The bottom panel shows that the F<2 g 
boron in-plane mode is very anharmonic. The inset indi- 
cates one of the E2 9 modes. Animations of these modes and 
their coupling with the electronic structure can be found at 



http: / /www.ncnr .nist.gov / staff / taner / mgb2 



The most interesting observation, though, is the totally 
opposite behavior of the E2 9 in-plane boron mode, as 
shown in Fig. 2b. The potential well is very shallow at 
small displacements and increases rapidly at large dis- 
placements, and can be fit very well to E(u) — A 2 u 2 + 
A4U 4 with a large ratio oiA^j A\ w 8. This indicates 
that this mode is unusually anharmonic and that we are 
in a non-perturbative regime. Hence the estimated har- 
monic energy, u>jj{Ei2g) = 60.3 meV will be lower than the 
actual energy. Using Hartree-Fock decoupling, one gets 
E(u) = (A 2 + 3A 4 < u 2 >)u 2 , where < u 2 >= /2Muj sch 
and uj sc h is the Self Consistent Harmonic (SCH) |E2| so- 



lution of ( 



sch 



[(A 2 +3A 4 < u 2 >)/M] 1 /2. This yields 
a better estimate of ui sc h ~ 70 meV, a 17% enhance- 
ment of the harmonic phonon energy u>H(E2 g ). We can 
calculate the exact energy levels of the potential by nu- 
merically solving the Schrodinger equation and obtain 
w{E2g) — 74.5 meV, a 25% enhancement of the harmonic 
value and also the Tc, which then matches with the peak 
in the GPDOS. 

Since the in-plane motions of the boron will change the 
boron orbital overlap, significant electron-phonon cou- 



pling can be expected for the planer a boron conduction 
band at the Fermi level and this plays an important role 
in the superconducting pairing. This can be most easily 
seen by comparing the band structure of the undistorted 
lattice to that of distorted one by a zone center phonon. 
For the harmonic phonons we did not see any significant 
changes near the Fermi level. However, for E2 9 modes, 
there is a significant splitting of the bands as well as shifts 
near the Fermi level as shown in Fig. 3, indicating strong 
electron-phonon coupling. 
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Fig. 3 Band structure of the undistorted (left) and dis- 
torted structures (right) by E2 9 phonons (its ~ 0.06A). Band 
structure of the other three frozen-phonon structures do not 
show any significant changes. 

To determine the electron-phonon (EP) coupling quan- 
titatively, we evaluated the Fermi-surface averaged defor- 
mation potential, 



A = ([Se(k)-S fJ ] 2 ), 



(1) 



for each zone-center frozen phonon p3] |. Here <5e(k) is the 
change in the one-particle energy with momentum k due 
to the frozen phonon, S/i the corresponding change in the 
chemical potential. <> denotes an average of k over the 
Fermi surface, which we have carried out numerically. 

For the harmonic Bi g , A2„, and Ei u phonons we cal- 
culated an insignificant coupling, and conclude that the 
electron-phonon coupling is negligible for these modes. 
For the F,2 g modes, Fig. 4a shows that the coupling A(u) 
is large and has both quadratic and quartic terms in 
frozen-phonon amplitude, u, indicating significant non- 
linear electron-phonon coupling. Neglecting non-linear 
cross terms, the electron-phonon coupling constant A for 
our anharmonic phonon is approximately given by [17] 



A = N(E F ) 



<n|Q|0> 
E„ — E 



< n\Q 2 \0 > 
E n — Eq 

(2) 



where N(Ep) is the total density of states at the Fermi 
energy and equal to 0.69 states/eV/cell. E n and \n > are 
the eigenvalues and eigenfunctions of the oscillator E2 ff 
in its adiabatic potential shown in Fig. 2b 



B' 2 and B\ 
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are the first and second-order EP coupling, respectively, 
and obtained from expansion of A(u). Q is the normal 
coordinate which is related to boron displacement u via 
a normalization constant. We calculate a total electron- 
phonon coupling A = 0.907(0.922) for M B =10 (11). Us- 
ing the McMillan expression for Tc p5| and taking a 
typical value for fx* = 0.15, we obtain a Tc of 39.4 K and 
38.6 K for Mb — 10 and Mb=11, respectively, yielding an 
isotope effect a — 0.21 |2(| (see Fig. 4b). These estimates 
are in excellent agreement with experiments. Within the 
above approximation p4j , the non-linear EP coupling in- 
creases Tc by about 10%. Our preliminary calculations 
of A including other non zone-center phonons indicates 
that the numbers will not change significantly. Therefore 
we conclude that MgB2 is a conventional electron-phonon 
superconductor in the strong-coupling regime. 



high Tc in MgB2. For hypothetical CaB2, the calcu- 
lated Abb is quite large (1.84 A) and we calculate very 
large anharmonicity, which may explain the instability of 
CaB2. It then appears likely that MgB2 is fortuitously 
just at the phase boundary. It will be interesting to see 
if it is possible to expand the <1bb in MgB 2 by Ca sub- 
stitution on the Mg site to increase Tc further. 

NOTE ADDED IN PROOF: After submission of our 
work, a recent Raman study |2£| has reported 77 meV 
for the Fi2 g mode with a very large width, confirming our 
predictions. 
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Fig. 4 (a) Calculated Fermi-surface averaged deformation 
potential A(u). The large quartic term indicates significant 
non-linear electron-phonon coupling, (b) T c versus boron 
mass, indicating a boron isotope effect a = 0.21, significantly 
reduced from harmonic BCS value of 0.5 due to giant anhar- 
monicity |2g| . Within our theory, we expect zero isotope effect 
for Mg, in good agreement with recent report of qm 9 =0.02 

It is interesting to note that there seems to be a close 
correlation between the anharmonicity of the E 2g in- 
plane modes and the B-B bond length (dss). For MgB 2 , 
dsB is 1.764 A, significantly stretched from its optimal 
value of 1.65 A in elemental planar boron, probably 
due to repulsive interactions between the Mg and B ions. 
This explains the unusual anharmonicity and observed 
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